home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_10 / a10_3.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  2.0 KB  |  92 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 10.3 (Crank-Nicholson Method for the Heat Equation).
  9. % Section    10.2,    Parabolic Equations, Page 517
  10. echo on; clc; format short; hold off; clear
  11.  
  12. % PARABOLIC EQUATIONS.
  13.  
  14. % Crank-Nicholson solution for the heat equation
  15.  
  16. %                          2
  17. %            u (x,t)   =  c  u  (x,t)
  18. %             t               xx
  19.  
  20. %  u(x,0) = f(x)   for  0 < x < a     and
  21.  
  22. %  u(0,t) = g1(t) and u(a,t) = g2(t)  for  0 ≤ t ≤ b
  23.  
  24. % A numerical approximation is computed over the rectangle
  25.  
  26. %         0 ≤ x ≤ a ,   0 ≤ t ≤ b.
  27.  
  28. pause    % Press any key to continue.
  29.  
  30. clc;
  31. %    Store f(x), g1(x), g2(x) in  f.m  g1.m  g2.m
  32. % function z = f(x)
  33. % z = sin(pi*x) + sin(3*pi*x);
  34.  
  35. % function z = g1(t)
  36. % z = 0;
  37.  
  38. % function z = g2(t)
  39. % z = 0;
  40.  
  41. delete f.m
  42. diary f.m; disp('function z = f(x)');...
  43.            disp('z = sin(pi*x) + sin(3*pi*x);');...
  44. diary off;
  45.  
  46. delete g1.m
  47. diary g1.m; disp('function z = g1(t)');...
  48.             disp('z = 0;');...
  49. diary off;
  50.  
  51. delete g2.m
  52. diary g2.m; disp('function z = g2(t)');...
  53.             disp('z = 0;');...
  54. diary off;
  55.  
  56. % Remark. f.m g1.m g2.m crnich.m trisys.m are used for Algorithm 10.3
  57. f(0); g1(0); g2(0);
  58. pause    % Press any key to continue.
  59.  
  60. clc;
  61. %    Place the endpoint of [0,a] in  a.
  62.  
  63. %    Place the endpoint of [0,b] in  b.
  64.  
  65. %    For the heat equation,  enter the constant  c.
  66.  
  67. %    Over [0,a] enter the number of grid points  n.
  68.  
  69. %    Over [0,b] enter the number of grid points  m.
  70.  
  71. a  =  1.0;
  72. b  =  0.1;
  73. c  =  1;
  74. n  =  11;
  75. m  =  11;
  76.  
  77. U = crnich('f','g1','g2',a,b,c,n,m);
  78.  
  79. pause    % Press any key to see the solution. 
  80.  
  81. Z = fliplr(U);
  82. mesh(Z);...
  83. Mx1 = 'The Crank-Nicholson solution to the heat equation.';...
  84. title(Mx1);...
  85. shg; pause % Press any key to continue.
  86.  
  87. W = U';
  88. points = W(:,2:n-1);
  89. clc,echo off,diary output,...
  90. disp(' '),disp(Mx1),disp(' '),disp(points),...
  91. diary off,echo on
  92.